Quantitative Molecular Orbital Energies within a GqWq Approximation 
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Using many-body perturbation theory within the GqWq approximation, we explore routes for 
computing the ionization potential (IP), electron affinity (EA), and fundamental gap of three gas- 
phase molecules - benzene, thiophene, and (1,4) diamino-benzene - and compare with experiments. 
We examine the dependence of the IP on the number of unoccupied states used to build the dielectric 
function and the self energy, as well as the dielectric function plane-wave cutoff. We find that 
with an effective completion strategy for approximating the unoccupied subspace, and a converged 
dielectric function kinetic energy cutoff, the computed IPs and EAs are in excellent quantitative 
agreement with available experiment (within 0.2 eV), indicating that a one-shot GqWo approach can 
be very accurate for calculating addition/removal energies of small organic molecules. Our results 
indicate that a sufficient dielectric function kinetic energy cutoff may be the limiting step for a wide 
application of Go Wo to larger organic systems. 
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INTRODUCTION 



Organic molecules and assemblies are of considerable interest for next-generation photovoltaics [IH3] and other 
energy conversion applications [JJ [5] . Their performance and utility hinges on understanding and control of their 
spectroscopic properties, such as ionization potentials (IPs) in gas-phase and solid-state environments, and orbital 
energy level alignment at interfaces. Density functional theory (DFT) is a widely used computational framework for 
studying structural and electronic properties of materials. However, Kohn-Sham frontier orbital energies and energy 
differences within common approximations to DFT, such as the local density approximation (LDA) and generalized 
gradient approximations (GGAs), are known to dramatically underestimate these quantities [BHS]. Recently, we have 
shown that accurate fundamental gaps for gas-phase and solid-state organic molecules [5] , and frontier orbital energies 
for an organic/metal interface [lOj ((1,4) diamino-benzene on Au(lll)) may be computed with many-body pertur- 
bation theory within the GW approximation 11.. For the latter, we found that the calculation must be adequately 
converged with respect to addition/removal energies of the isolated components, i.e. molecule and substrate. In this 
article, building on prior work [TUJ IT2Tf2T)] . we explore the extent to which we may obtain accurate IPs and electron 
affinities (EAs) of gas-phase molecules using a GqWq approximation. 

While there are numerous studies benchmarking the GW approximation against transport gaps in bulk inorganic 
solids [111 I21f[2"4l . similar works for isolated molecular systems are less common, and while all works exhibit marked 
improvement over standard DFT approaches, there is some quantitative disagreement (see e.g., [TUJ [12 19. 25]). For 
example, for gas-phase molecules, using an atom-centered basis set, it has been found that self-consistency in either 
the GW eigenvalues [12 -14 or in both the eigenvalues and eigenvectors [T5] is essential for obtaining good agreement 
of computed molecular IP and EA with experiment. On the other hand, with a planewave basis set, it has been 
demonstrated that a more systematic representation of the dielectric matrix and Coulomb- hole (CH) term, £ch, 
brings the GoHVpredicted IP and EA in closer agreement with experiment [TUJ [TUl 02] ■ Beyond the differences in 
their basis sets, these studies have differed in their representation of the dielectric matrix, the presence of a truncation 
scheme for the Coulomb interaction, and their approach for handling the empty states necessary to converge Sqh • As 
a consequence, the accuracy of different GW approaches for the IP of gas-phase molecules remains an open question. 

Here, we compute the GqWq IP, EA, and fundamental gap (IP - EA) of three gas-phase molecules benzene (BEN), 
(1,4) diamino-benzene (BDA), and thiophene (TP), as shown in Fig. [I] and compare the computed IP and EA with 
measurements |26H29j . We examine the dependence of the IP and fundamental gap on the number of unoccupied 
states used to build the dielectric function and the self energy, as well as the dielectric function G-space cutoff. We 
find that as our calculations approach convergence, the computed IPs and EAs are in excellent quantitative agreement 
with experiment (within 0.2 eV), indicating that GqWq can be very accurate for calculating addition/removal energies 
of small organic molecules. 
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METHODS 

Our GW calculations are performed using the BerkeleyGW [30 package, following an established GqWo ap- 
proach [TT] . The self-energy, E = iGW, is computed as a first order correction to the Kohn-Sham DFT Hamiltonian. 
The quasiparticle states are taken from DFT within the GGA of Perdew, Burke, and Ernzerhof (PBE) [5T] and are 
expanded in a planewave basis set. The cutoff for the planewave expansion is 80 Ry for BEN and TP and 60 Ry for 
BDA, and is determined such that the DFT total energy is converged to < 1 meV/atom. The molecular geometry is 
optimized such that forces are less than .04 eV/A. Norm conserving pseudopotentials are used, with 1, 4, 5, and 6 
electrons explicitly treated as valence for H, C, N, and S, respectively. 

Since periodic boundary conditions are imposed in our planewave DFT and subsequent GW calculations, the 
molecules are placed in a large supercell, chosen to be twice the size necessary to contain > 99% of their charge 
density. The supercell dimensions are 14 x 8 x 15 A 3 for BEN, 14 x 9 x 14 A 3 for TP, and 15 x 15 x 15 A 3 for BDA. 
In constructing the dielectric matrix and the self-energy at the Go Wo step, the Coulomb potential is truncated at half 
of the unit cell length in order to avoid spurious interactions between periodic images. The electrostatic potential 
at the surface of the supercell is computed at the DFT level and its average subtracted from the GW eigenvalues to 
obtain absolute energies and therefore, IPs and EAs. 

The static inverse dielectric function (cg 1 g'( c i)) * s expanded in planewaves (with wavector G), and a cutoff (e^* = 
|q+G| 2 /2), where q is a wavevector. e^Q, (q) is constructed as a sum over unoccupied states [32 , which is truncated 
at a finite number of states, N c , with energy E(N e ). The dielectric function is extended to finite frequency with the 
generalized plasmon-pole (GPP) model of Hybertsen and Louie [55] . 

For the purposes of analysis, we define the self-energy operator as a sum of Fock exchange, screened exchange, and 
the CH terms [llj . The screened exchange term, Esx, requires an explicit sum over just occupied states; however, it 
is implicitly dependent on N e through the dielectric function. The CH term, Sch, involves a sum which in principle 
must span the full unoccupied subspace, but in practice is also truncated at finite number of unoccupied states N c , 
with corresponding energy E(A C ). For simplicity, we set N e equal to N c , subtract the matrix elements of the Fock 
operator, Ex (which is independent of N c and e), from Esx, and study the convergence behavior of Esx-x and Ech 
terms with respect to N c and £q ! . 

RESULTS 
Convergence of the dielectric matrix 

Fig. [l] summarizes our calculated IPs for BEN, TP, and BDA as a function of two parameters, £q* and N c . The 
IP is defined here such that a positive value indicates a bound electron. The IP increases significantly, by about 0.5 
eV, as either parameter is increased (taken towards convergence) for all three molecules; in contrast, the fundamental 
gap, IP-EA, converges rapidly to within 0.1 eV for E(N C ) > 2 Ry and £q* > 4 Ry. 
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FIG. 1. The GqWq predicted IP and fundamental gap as a function of e 1 cutoff and N c . 



As noted previously [33], the interdependence of N c and Cq* can lead to a "false convergence" of the IP with 
respect to the dielectric cutoff at small fixed N c . For all three molecules at E(N C ) = 2 Ry and £q* = 4 Ry, the IP 
is apparently converged to within 0.1 eV; however, if E(N C ) is increased to 6 Ry, the IP varies by 0.3 — 0.4 eV as 
6q* is raised from 4 to 24 Ry. For all three molecules studied, this "false convergence" subsides for E(N C ) ~ 6 Ry 
above the vacuum level (corresponding to N c ~ 3000 for BEN and TP, and ~ 5000 for BDA within our supercells); 
the computed IP is unaffected by an increase of £q 1 for values greater than 12 Ry for E(N C ) > 6 Ry. However, the 
IP is still quite sensitive to N c , as we will discuss further below. 

For fixed N c , both Egx-x and Ech also appear converged (to within 0.1 eV) for £q ! > 12 Ry, as shown in Fig. [2] 
for the highest occupied molecular orbitals (HOMOs) of BEN, TP, and BDA. Interestingly, for low N c the variation 
of Sgx-x and Ech with £q * ranges from 2 — 100 times larger than the corresponding variation of the IP. Thus, Esx-x 
and Ech are evidently less prone to "false convergence" at low N c than the IP. Since both Esx-x and Ech depend on 
e , but with opposite sign 32!, their sum (which determines the IP) is less sensitive to an underconverged dielectric 
function. 

While Fig. [l] and Fig. [2] suggest that an £q' > 12 Ry is sufficient for a precision of 0.1 eV or better in the IP for 
fixed N c , they also highlight the fact that the self-energy corrections are more sensitive to N c when the high energy 
Fourier components of are well-described. Fig. [l] shows a variation in IP of > 1 eV as E(AT C ) grows from 2 Ry 
to greater than 6 Ry. Fig. [2] indicates that the Ech term is responsible for this variation, as E$x appears converged 
within 0.2 eV for a dielectric matrix described with E(N C ) > 6 Ry. This implies that for the molecules and supercells 
under study, for £q' > 12 Ry and E{N C ) > 6 Ry, the only remaining convergence issue in the calculation is the sum 
over the unoccupied subspace. We now discuss the different strategies for converging this sum. 
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FIG. 2. For BEN, TP, and BDA HOMOs: E sx _ x (solid lines) and E CH (dashed lines) as a function N c and eg 1 * for E(N C ) = 
25, 50, 80, and 105 eV. The legend follows Fig. 1. 
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FIG. 3. a) Ech as a function of number of bands for the BDA HOMO for both static COHSEX and G W . The static 
COHSEX result for N c — > oo is indicated with a horizontal dotted line, b) The GoWo IP with the CH term extrapolated to 
infinite N c using fitting techniques and the static remainder approach. 



Convergence of the Coulomb-hole term of the self-energy 

The slow convergence of the Ech term, for a converged value of £q * , with respect to N c can be seen in Fig. [3ja. for the 
BDA HOMO. E C H varies by more than 2 eV for N c e [500; 5000] and shows a finite slope of 10~ 4 eV/N c at N c = 5000. 
Moreover, this same slow convergence behavior can be seen with a static CH and screened exchange method (static 
COHSEX) for which a full evaluation (shown as dashed line) does not require a sum of empty states [34] . Comparison 
of our dynamic and static calculations suggests that the N c dependence of Ech comes from both static and dynamical 
correlation terms. The static COHSEX CH term is still 0.2 eV away from the exact solution at N c — 5000, and has 
a different slope than the full dynamical Ech- 

The slow convergence of Ech with N c has been addressed with different strategies in prior work [T§1 l3"5Tl4"0] . Here, 
we examine three different approaches for extrapolating the CH term to infinite N c and examine their consequence for 
the IP: i) fitting Ech (-^c) f° r a given orbital with an analytical form, and calculating its limit when N c — > oo (see, e.g. 
[H]); it) fitting the dynamical Ech (N c ) to a functional form determined from the corresponding static COHSEX 
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term |10J ; and Hi) approximating the correction to the dynamical CH term based on completing the unoccupied 
subspace within the static COHSEX approximation, i.e. the static remainder (SR) approach [36 . 

Kang and Hybertsen applied a fitting scheme to £ch to obtain the valence band maximum of Ti02 and found a 
0.2 eV range in predicted values for two different functional forms for the fit [H]. We take a similar approach and 
consider the following four functional forms for the dynamical Sch (-^c) : 



£ CH (Ay ~a + /3N^, (1) 

£ CH (Ay ~a + pN~\ (2) 

Ecn(N c )~a + PN-z, (3) 

X CR (N c )~a + (3e--, (4) 



where a, /?, and 7 are fitting parameters. In practice, we find that good fits (P value < 0.005) can be consistently 
obtained using any of these forms. 

We also fit the partial sum Sch(Ac) computed within static COHSEX such that a is the numerically exact closed 
form value of the static CH (Egg tic (oo)). More precisely, the static CH term, Sgg tic (N c ), is fit to Eq. [l] with (3 and 
7 as fitting parameters. The dynamical Y>cn(N c ) is then fit to Eq. [Ij with 7 fixed and a and (3 as fitting parameters. 
Here, we are assuming that the same functional form describing the static £ch also describes the dynamical case. 

Lastly, we apply the SR correction defined in Ref. |36j where 

S CH (N c -> 00) ~ E CH (N c ) + l - [Eg| tic (oo) - £ s c ^ ic (N c )} . (5) 

In Fig. |3Jd) , we report the computed IPs of BDA using all five extrapolation techniques described above. Because 
we are far from convergence in N c , the fitting procedure (i) is much less favorable than found by Ref. [41j both by its 
error with respect to experiments and its range of uncertainty: the assigned functional form can produce predicted 
IPs ranging from 5.8 to 7.2 eV. More importantly, the computed IP is very sensitive to the number of bands initially 
used. The best fit to the static COHSEX result for Sch, (H), results in IPs that monotonically increase with the 
number of bands used in the fit, and appears to be converging towards the SR result. 

The SR method gives the best results, with predicted IP values within 0.1 eV for N c g [500; 5000]. The results of 
the SR method are particularly remarkable in the sense that when using this procedure, less unoccupied states are 
needed to converge the CH term than the dielectric matrix (respectively 500 and 5000 for BDA). 

Comparison with experiment 

Table [T] shows the GqWq IP and EA for BEN, TP, and BDA, along with experimental values. For all molecules, 
we use E(N C ) = 6 Ry and £q* = 24 Ry, and Sqh is extrapolated to infinite number of bands via SR [3^. Our 
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TABLE I. Go Wo IP for BEN, TP, and BDA in eV. The calculations are performed with eg" or 24 Ry, with E(N C ) fixed at 6 
Ry, and the static remainder correction applied. 



Molecule 


BEN 


TP 


BDA 


IP Theory 


9.4 


9.0 


7.3 


IP experiment 


9.24 [26] 


8.86 [26] 


7.34 [13[28] 


EA Theory 


-0.92 


-0.94 


-0.90 


EA experiment 


-1.1 [29] 


-m 





Go Wo results are in excellent agreement with experiment, within 0.2 eV for IP of all three molecules and the EA of 
BEN. Our predictions agree well with previous planewave-based Go Wo studies [17TH9] for BEN, but differ somewhat 
quantitatively with with other Go Wo results obtained using localized basis sets for TP [T3] and BDA [2D] . 

CONCLUSIONS 

With use of unoccupied states that span ~ 6 Ry in energy, an e^' greater than or equal to 12 Ry, and the static 
remainder approach to correct for the finite number of empty states in Sch, we obtain converged values for the 
GoWo-calculated IP and EA of three organic molecules in the gas-phase. The predicted IPs and EAs agree to within 
0.2 eV with available experiment. Our results indicate that GoWo provide quantitatively accurate addition/removal 
energies for small organic molecules. We find that a limiting step to these calculations is the large £q* required for 
convergence. Thus, extrapolation techniques for £q* will be increasingly valuable for describing larger systems, such 
as metal/organic molecule interfaces. 
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